A Statistical Social Network Model for Consumption Data in Food Webs

We adapt existing statistical modelling techniques for social networks to study consumption data observed in food webs. These data describe the feeding among organisms grouped into nodes that form the food web. Model complexity arises due to the extensive amount of zeros in the data, as each node in the web is predator / prey to only a small number of other nodes. Many of the zeros are regarded as structural (non-random) in the context of feeding behaviour. The presence of “bottom prey” and “top predator” species (those who never consume and those who are never consumed, with probability 1) creates additional complexity to the modelling framework. We develop a special statistical social network model to account for such network features. The model is applied to a well-studied food web, for which the population size of seals is of concern to various commercial fisheries. 1Corresponding author. CSIRO Mathematics, Informatics and Statistics (CMIS), GPO Box 664, Canberra, ACT 2601, Australia. E-mail: [email protected] 2Department of Mathematical Sciences, University of Nevada, Las Vegas, NV 89154, USA. E-mail: [email protected] ar X iv :1 00 6. 44 32 v2 [ st at .M E ] 2 4 A ug 2 01 0 1 Why social networks? A food web is a network of organisms. When the relationship among them is of interest, organisms are typically aggregated at various resolutions to form nodes. For example, one node may consist of various squid species, while another may consist of the single species Argyrosomus hololepidotus (kob). In the trophic context, we are interested in the feeding relations among the nodes in the food web. Consider the pair of nodes (i, j). The within-pair predatory relation can be depicted as one of the following: i j i→ j i← j i↔ j (1) where, conventionally, any link / arrow points from prey to predator. From left to right, the depictions in (1) respectively represent no predation between i and j, predation of i by j but not vice versa, predation of j by i but not vice versa, and mutual predation between i and j. To represent (1) in a quantitative framework, the absence of a link is referred to as a zero link, so that “i j” consists of two zero links, each of “i→ j” and “i← j” consists of a zero and a non-zero link, and “i ↔ j” consists of two non-zero links. Since these links are directed, each pair (i, j) yields two directed links: from i to j, and from j to i. Extending this to all n nodes in the food web, we have a network that consists of 2×(n-choose-2) = n(n− 1) pairwise or dyadic directed links. Research on relational patterns in a network of nodes arises in many practical settings, most commonly in the social sciences, e.g. friendship, international trade. In such contexts, networks are referred to as social networks. We apply the same nomenclature to studying relational patterns in a food web with respect to predation. Various quantitative social network analysis (SNA) techniques have been developed to understand network relational patterns (e.g. Mucha et al., 2010; Wasserman & Faust, 1994) and adopted in food web research (e.g. Dambacher et al., in press; Krause et al., 2003); these SNA methods are based largely on the mathematical notion of equivalence class for defining congruence among individual elements in a given set according to certain criteria (see, e.g. Düntsch and Gediga, 2000). The objective of this type of SNA is to seek optimal partitions of the network into compartments of nodes subject to the given criteria. For example, compartments identified in a food web may correspond to trophic levels. Recently, statistical regression methodologies have been developed to express network links as the random response of within-node and internode characteristics (Westveld & Hoff, in press; Hoff, 2005; Gill & Swartz, 2001). The resulting conclusions about network features are purely empirical, based entirely on observed network attributes without the use of network dynamics theory or subject-matter theory. Chiu & Westveld (2010) demonstrate that, in the context of food webs, statistical SNA techniques can provide an alternative perspective of trophic relational patterns according to feeding activity and preference. These authors apply the statistical SNA modelling framework by Ward et al. (2007) to regress the presence-absence of pairwise predation, in a mixed-effects logistic model, on the dyadic node characteristic of phylogenetic similarity; eight food webs are analyzed this way. The basis of a statistical model for SNA is a two-way analysis-of-(co)variance (ANO(CO)VA). To see this, consider the mixed model in Chiu & Westveld (2010) with one optional covariate: log pij 1− pij = β0 + β1xij + si + rj + u ′ ivj + εij , i 6= j (2) 2 where pij is the probability that j predates on i, xij = xji is the phylogenetic similarity between i and j, si is the ith random sender (prey) effect, rj is the jth random receiver (predator) effect, uivj is the random interaction between i and j expressed as an inner product of k-dimensional vectors ui and vj , and εij is the remaining random component unattributable to xij, si, rj,ui, or vj . All of si, rj,ui,vj, and εij are mean-zero Gaussian random errors. (Note that expressing the interaction term as the inner product of latent vectors ui and vj is due to Hoff et al., 2002. The latent uand v-spaces are abstract entities; their dimension k can be regarded as a model parameter to be estimated from the data. See Section 1.1 for an interpretation of the latent spaces.) Complex network dependence not addressed by the ANO(CO)VA equation is modelled through Cov(si, ri) = Σ = [ σ s ρsrσsσr ρsrσsσr σ 2 r ] for all i = 1, . . . , n , (3) Cov(εij, εji) = Ω = σ [ 1 ρ ρ 1 ] for all i 6= j . (4) Equation (3) stipulates that sender and receiver effects due to the same node are possibly correlated (ρsr), with potentially distinct amounts of uncertainty (σs and σr); this within-node structure is constant across nodes. Equation (4) allows for potential reciprocity between i and j through ρ (constant across all (i, j)-pairs), with the typical assumption of homogeneous random errors through σ. There are various advantages to the use of a statistical SNA framework when studying food webs: 1. Randomness or uncertainty is naturally inherent in the network links. For example, having observed “i j” (no predation) does not necessarily preclude the future occurrence of “i→ j.” This type of randomness is readily acknowledged and modelled through the statistical framework. 2. Factors that may influence or be associated with the network links can be explicitly expressed as covariates in the statistical model. This facilitates the immediate understanding of what makes a given node the predator or the prey, and avoids post-SNA “detective work” that may be necessary when using non-statistical SNA methods, which provide information about who tends to eat whom but not why. Depending on the context, covariates in a statistical model often provide insight into why. 3. Statistical inference for the random effects si, rj,ui, and vj provides insight into trophic patterns (e.g. clustering of nodes) from three different perspectives simultaneously: (i) activity level as prey and as predator, (ii) preference of being consumed, and (iii) preference of consuming; see Section 1.1 for the rationale. A single SNA that considers all three perspectives provides a more comprehensive understanding of food web structure. Assessing uncertainty in the trophic patterns is part of the one-step inference. With equivalence class methods, typically one analysis is produced for each pre-defined set of criteria corresponding to a single perspective of congruence, and uncertainty assessment is applied separately by randomly permuting network links then re-analyzing each permutation under the given criteria (Krause et al., 2003). 3 4. As part of the one-step inference through a model fit, statistical inference for network dependency parameters ρsr and ρ provides further insight into trophic features that underlie the link patterns (see Section 1.1). 5. The regression framework of a statistical SNA readily allows (i) prediction of network links under different scenarios, and (ii) assessment of uncertainty in the predictions. 6. Since a statistical SNA is based solely on empirical information, results in certain contexts may be used to validate projections made based on deterministic mathematical models for trophic relations. 1.1 Interpreting the statistical social network model for food webs Although observations on feeding behaviour may be observed at the organism or species level, trophic food web data are typically recorded for nodes and presented in a diet matrix, as follows: Node as Node as Predator Prey 1 2 . . . n 1 y11 y12 . . . y1n 2 y21 y22 . . . y2n .. .. .. .. .. n yn1 yn2 . . . ynn For presence-absence data, yij = 1 if j is observed to predate on i, and yij = 0 if this predation is not observed. For weighted data, yij is the weight, or magnitude, of the consumption of i by j. For example, yij can be the total biomass of i consumed by j. Whether presence-absence or magnitude is considered, yii corresponds to cannibalism of i, and is not modelled by existing statistical SNA methods. This is because “i → i” is often undefined in a typical social network context such as friendship or trade. Although cannibalism is not uncommon in nature, we regard it as a secondary or nuisance feature of the food web when internodal trophic relations are considered empirically. The reason is as follows. Based purely on which node is observed to predate on which other node, we wish to understand primarily the relative influence of nodes on each other, instead of a given node’s self influence. Furthermore, a model such as (2) depicts a temporal snapshot of the food web, and does not address food web dynamics in which cannibalistic behaviour may directly influence other nodes due to the dynamical nature of the system. An empirical analog of a dynamic network model is the longitudinal social network (LSN) model of Westveld & Hoff (in press). While potentially valuable in principle for food web research, longitudinal statistical SNA often may be infeasible in this context due to the lack of repeated field observations of the same food web over time. Consequently, in this report we focus on the non-temporal empirical aspect of food webs while ignoring cannibalism. Despite the limitations of the statistical SNA framework, it can provide valuable insight into the food web structure through its greatly interpretable model parameters. To see this, consider (2)–(4). Given Node i, the bivariate random effect [si, ri] describes the level of feeding activity of the node, after adjusting for covariate effects. Feeding activity for i is related to its in-degree 4 ( ∑ j yji = total activity as predator) and out-degree ( ∑ i yij = total activity as prey). Inand outdegree unexplained by covariate(s) are captured respectively by si and ri. A graph of the estimated [si, ri] ′ vectors (Graph SR in Chiu & Westveld, 2010) displays the positions of nodes according to their feeding activity; any identifiable cluster may be considered a “trophic level” from the perspective of feeding activity as prey and as predator. Estimates of the k-dimensional vectors ui and vj may also be graphed (Graphs U and V in Chiu & Westveld, 2010). We take k = 2 to reduce model complexity and allow easy visualization of the vectors. Alternatively, k can be determined via optimization criteria, then projected onto R for graphical display (Ward et al., 2007; Hoff, 2005). In either case, the latent k-dimensional u-space corresponds to preference of being consumed. As Chiu & Westveld (2010) explain, if ui and uj are neighbours in the u-space, then the sending behaviour of i to ` — after accounting for sending activity — is similar to that of j to `, for all nodes `. In a food web, this phenomenon translates to i and j being similarly preferred as prey. The same interpretation applies to vi and vj being neighbours in the v-space, except for their similarity in receiving, or preference for prey. Thus, clustering in the u-space suggests trophic levels with respect to preference of being prey, and clustering in the v-space suggests trophic levels with respect to preference of being predator. Chiu & Westveld (2010) demonstrate that trophic clusters identified in the three graphs can differ substantially depending on the perspective from which trophic relations are viewed. For example, nodes which show similarity in feeding activity ([si, ri] vectors being closely clustered) may differ drastically in their feeding preference (vi vectors belonging to separate clusters in the v-space); this was seen in various food webs analyzed by Chiu & Westveld (2010). Finally, statistical SNA not only provides information on trophic clustering from various perspectives, but also on the dependency among nodes beyond the patterns of internodal links. This extra insight is achieved through the statistical inference for ρsr and ρ. Consider the sign of ρsr, which is reflected by the trend of the graphical display of the estimated [si, ri] vectors. A positive trend suggests that active predators tend to be active as prey, and a negative trend suggests that active predators tend to be inactive as prey. Insight into the reciprocity of predation, or the tendency of predator-prey role reversal, is available through the inference for ρ. A positive ρ indicates that the predation of i by j is positively associated with the predation of j by i, and hence, the tendency for predator-prey role reversal between i and j. Conversely, a negative ρ indicates that this role reversal is unlikely. Note that statistical social network models typically involve Bayesian inference, as mixedeffects models with complex dependence structure can be naturally constructed in a Bayesian hierarchical framework. Posterior inference can be made via Markov chain Monte Carlo (MCMC). However, posterior inference for ui and vj via MCMC requires careful handling of the MCMC samples; see Chiu & Westveld (2010). 2 Modelling consumption data The discussion on the method of statistical SNA and its merits apply generally to the cases of yij being presence-absence or weighted data. However, weighted food web data may pose a challenge, due to the high incidence of yij = 0. For example, the eight food webs analyzed by Chiu & 5 Westveld (2010) each consists of between 76% and 98% zeros. Direct application of existing statistical SNA techniques to such weighted data would require a special distributional assumption for yij to account for its extreme point mass at 0 and its continuous distribution away from 0 (Figure 1). For this, a mixture distribution may be appropriate, at the expense of model complexity and computational burden. In this report, we propose an alternative approach that does not require the same level of model complexity, and is reasonably straightforward to implement. Figure 1: Distribution of consumption data for the Benguela food web. Population Consumption (20th Root) F re qu en cy 0.0 0.5 1.0 1.5 0 10 0 20 0 30 0 40 0 50 0 60 0 Consider the nature of zeros in food web data. It is unlike the social context of, say, friendship in which typically (i) non-zero links between two nodes are common, and (ii) the randomness inherent in the links is substantial enough that under different scenarios, zero links can plausibly become non-zero and vice versa. In contrast, zero links often dominate a food web; given such a link, biological theory can easily identify the nature of this 0, such as the case of herbivorous grazers almost surely not consuming other animal species under any realistic scenario. Among the vast number of zeros, those regarded as truly random are typically no more than a handful, if any. For this reason, here we regard all observed zeros as structural zeros, and remove them from consideration when constructing the statistical social network model. Note that for presence-absence food web data, including the zeros in the logistic model (as is done by Chiu & Westveld, 2010) implies a different interpretation than what we propose in this report. Specifically, (2) regards the entire food web (i.e. the set of all n(n − 1) directed links) as one random entity, whereas in our current context, the randomness in each directed link is being modelled. As discussed above, it is not straightforward to consider the randomness of the entire weighted food web, for which the 6 distribution of yij is highly non-standard. We refer to the weighted food web excluding zero links as the “reduced weighted food web,” or “reduced web” for short. To model this reduced web, some notation is necessary. Let S∗ = {(i, j) : j > i, i, j = 1, . . . , n} , S0 = {(i, j) ∈ S∗ : yij, yji = 0} , S = S∗ \ S0 , S1 = {(i, j) ∈ S : yij > 0, yji > 0} , S2 = {(i, j) ∈ S : yij > 0, yji = 0} , S3 = {(i, j) ∈ S : yij = 0, yji > 0} . Note that S∗ is the set of all n-choose-2 nodes in the food web, and S0,S1,S2,S3 are disjoint. We need not consider the set S0 of unlinked pairs, since by definition of the reduced web, such pairs are not considered. Thus, ∪k=1Sk = S. Also note that S1 consists of all pairs who show mutual predation, S2 consists of send-only pairs (i.e. j predates on i but not vice versa), and S3 consists of receive-only pairs (i.e. i predates on j but not vice versa). We additionally let I = {1, 2, . . . , n} , I1 = { i ∈ I : n ∑ j=1, j 6=i yij > 0, n ∑ j=1, j 6=i yji = 0 } , I2 = { i ∈ I : n ∑ j=1, j 6=i yij > 0, n ∑ j=1, j 6=i yji > 0 } , I3 = { i ∈ I : n ∑ j=1, j 6=i yij = 0, n ∑ j=1, j 6=i yji > 0 } . The set I consists of all n nodes of the food web, which is broken down into three disjoint sets I1, I2, and I3, where ∪k=1Ik = I. The set I1 consists of all basal nodes, or bottom prey who never consume but are predated upon by at least one other node in the web; in contrast, I3 consists of top predators who only consume but are never predated upon; I2 consists of the remaining “middle” nodes. For a given node, the standard statistical social network model simultaneously considers links into and out of the node through (3) and (4). For the reduced web, however, (3) only applies to nodes in I2, each of which plays the role of prey as well as that of predator. In contrast, nodes in I1 and I3 are linked to other nodes in one direction only, so that (4) is degenerate and reduces to si being independent and identically distributed (iid) as N(0, σ s ) for i ∈ I1 and ri being iid N(0, σ r ) for i ∈ I3. Similarly, (4) only applies to mutual predators (i, j) ∈ S1. For S2 or S3, one of the two links is missing, so that (4) is degenerate and reduces to εij being iid N(0, σ). Under this reduced framework, technically si is undefined for i ∈ I3 and rj is undefined for i ∈ I1. However, to facilitate the visualization of feeding activity, we arbitrarily define si ≡ −4σ s for all i ∈ I3 , ri ≡ −4σ r for all i ∈ I1 (5) 7 so that given the variability of overall feeding activity, the “random effects” in (5) are in fact constant, and are appropriately located in the far left tail of the distribution of si and rj . Then, any [si, ri] ′ for i ∈ I may be displayed on the sr-plane. On the other hand, leaving ui,vi for i / ∈ I2 undefined would not hinder the visualization of feeding preference, since we consider the uand v-spaces separately (as opposed to the cross product of the sand r-spaces). Thus, we can consider the distribution of ui for all i ∈ I1 ∪ I2 in the u-space, and that of vi for all i ∈ I2 ∪ I3 in the v-space. As do Chiu & Westveld (2010), we also take k = 2 to minimize model complexity, and take Var(uiq) = σ uq, and Var(viq) = σ 2 vq for all i and q = 1, 2, where ui = [ui1, ui2] and vi = [vi1, vi2]. 3 Model implementation for the Benguela food web A well-studied food web is that of the Benguela ecosystem, originally discussed by Yodzis (1998) but further studied by, e.g. Chiu & Westveld (2010) and Dunne et al. (2004). 3.1 Data The diet matrix in Yodzis (1998) consists of diet percentages for 29 nodes (Table 1), and is accompanied by several relevant variables: adult individual body mass, annual harvest, carrying capacity, ingestion factor, and population biomass. The (i, j)th cell of the diet matrix represents the percentage of j’s diet through consuming i. Thus, each column of the diet matrix necessarily sums to 100 if all organisms in the actual food web are represented in the diet matrix. However, this appears not to be the case, as column sums range from 90 to 131.5 for non-basal nodes. To derive weighted data that corresponds to consumption volume, first we scale the diet percentages according to the column sums so that each scaled column for a non-basal node sums to 100. The columns for the two basal nodes necessarily consist of all 0s and remain unaltered. Next, we use the scaled diet proportions to define two different measures of consumption volume: (i, j)th population consumption = (j’s population biomass)× (j’s ingestion factor)× (j’s diet proportion due to i) (i, j)th per-adult consumption = (j’s adult individual biomass)× (j’s ingestion factor)× (j’s diet proportion due to i) where j’s ingestion factor is the “fraction of physiologically maximal ingestion” (Yodzis, 1998) for j. Both definitions refer to the biomass of i consumed by j. We take yij = ((i, j)th population consumption) 1 20 (6) or yij = ((i, j)th per-adult consumption) 1 10 . (7) 8 As can be seen in Figure 1, the 20th-root transformation in (6) results in a reasonably Gaussian distribution after the removal of 0 links. The same is true for the 10th-root transformation in (7) (not shown). The objective is to express either definition of yij as the response of covariates in the reduced statistical SNA framework. Although these transformations appear to be drastic and not easily interpretable, having approximate normality eliminates extra model complexity that would have been needed to address non-normality, given an already complex framework for modelling network dependence. 3.2 Bayesian hierarchical model For these data, only 196 out of the 29×28 = 812 pairwise links are non-zero. The sets S3 = I3 = ∅ as there are no receive-only nodes. The cardinality of S1,S2, I1, and I2 are, respectively, 5 (see the illustration on the cover of this report), 186, 2, and 27. Thus, the statistical social network model comprises the following set of distributional statements: [ yij yji ]∣∣∣∣μij, μji, si, sj, ri, rj,ui,uj,vi,vj,Ω ∼ BVN( [ μij μji ] + [ si ri ] + [ rj sj ] + [ uivj ujvi ] ,Ω) for all (i, j) ∈ S1 , (8) yij|μij, si, rj,ui,vj, σ ∼ N(μij + si + rj + uivj, σ) for all (i, j) ∈ S2 , μij = x ′ ijβ for all (i, j) ∈ S , si|σs ∼ N(0, σ s) , ri|σr ≡ −4σr for all i ∈ I1 , [ si ri ]∣∣∣∣Σ ∼ BVN(0,Σ) for all i ∈ I2 , (9) uiq|σuq ∼ N(0, σ uq) for all i ∈ I and q = 1, 2 , viq|σvq ∼ N(0, σ vq) for all i ∈ I2 and q = 1, 2 , where xij is the vector of covariate values for pair (i, j), and β is the corresponding regression coefficient. The choice of xij is discussed in Section 3.3. For Bayesian inference of complex models, the mixing of MCMC draws is often a practical concern. To mitigate mixing difficulties, first note that (9) implies si ∼ N(0, σ s) , ri|si ∼ N(λsi, φ) , σ r = φ + λσ s , ρsr = λ σs σr . (10) To see this, rewrite the last expression of (10) as λ = ρsrσr/σs. Hence, the 2nd expression of (10) implies φ = (1− ρsr)σ r . These expressions for λ and φ imply ρsr = λ 2σ 2 s σ2 r = σ r − φ σ2 r which in turn implies the third expression in (10). The use of (10) avoids generating MCMC samples of Σ from a matrix distribution such as Wishart, which caused major mixing difficulties 9 in our case. Mixing issues associated with Wishart priors are also discussed in The BUGS Project frequently asked questions. We consider reasonably diffuse proper priors in the Bayesian hierarchy. The diffuseness reflects our lack of prior knowledge of parameter values. Let Γ(a, b) denote the Gamma distribution parametrized in such a way that small values of a and b lead to diffuseness. We take ρ = e − 1 e2z + 1 , z ∼ N(0, 0.82) , (11) λ, β` ∼ N(0, a−1) for all ` = 0, 1, . . . , L , (12) σ−2, σ−2 s , φ −2, σ−2 uq , σ −2 vq ∼ Γ(a, a) for all q = 1, 2 , where L is the number of covariates, and a < 0.01 (small a leads to diffuseness). The actual choice of a varied in our analyses and was found to have virtually no influence on the results. Expression (11) employs the Fisher transformation to avoid taking ρ ∼ U[−1, 1] so as to improve MCMC mixing; the choice of distribution for z corresponds very closely to ρ ∼ U[−1, 1]. 3.3 Covariates From the variables in Yodzis (1998) that accompany the matrix of diet percentages, we consider ti1 = sender adult individual biomass, ti2 = annual harvest of sender, tj1 = receiver adult individual biomass, tj2 = annual harvest of receiver. Note that ts are sender-specific covariates, and ts are receiver-specific. In addition, we create a pair-specific covariate tpij = taxonomic distance by comparing the taxonomic classification of i and j according to their domain, kingdom, phylum, class, order, family, genus, and species. This is the distance counterpart of the conservative phylogenetic similarity measure zij of Chiu & Westveld (2010). We do not consider carrying capacity as a covariate, as it is highly related to the notion of ingestion factor, which is used to define the response variable. Altogether, we have five covariates, all of which are non-negative. A log-transformation is applied to all five, which substantially reduces skewness of the distributions. Covariates that have 0 values are shifted up by the value 1 before taking log, and thus non-negativity is preserved. Next, each log(t) is centered to reduce dependence among the corresponding regression coefficients in the Bayesian inference. We assume this dependence to be negligible, as is reflected by